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Efficient techniques to navigate networks with local information are fundamental to sample large- 
scale online social systems and to retrieve resources in peer-to-peer systems. Biased random walks, 
i.e. walks whose motion is biased on properties of neighbouring nodes, have been largely exploited 
to design smart local strategies to explore a network, for instance by constructing maximally mixing 
trajectories or by allowing an almost uniform sampling of the nodes. Here we introduce and study 
biased random walks on multiplex networks, graphs where the nodes are related through different 
types of links organised in distinct and interacting layers, and we provide analytical solutions for 
their long-time properties, including the stationary occupation probability distribution and the 
entropy rate. We focus on degree-biased random walks and distinguish between two classes of walks, 
namely those whose transition probability depends on a number of parameters which is extensive 
in the number of layers, and those whose motion depends on intrinsically multiplex properties 
of the neighbouring nodes. We analyse the effect of the structure of the multiplex network on 
the steady-state behaviour of the walkers, and we find that heterogeneous degree distributions as 
well as the presence of inter-layer degree correlations and edge overlap determine the extent to 
which a multiplex can be efficiently explored by a biased walk. Finally we show that, in real-world 
multiplex transportation networks, the trade-off between efficient navigation and resilience to link 
failure has resulted into systems whose diffusion properties are qualitatively different from those of 
appropriately randomised multiplex graphs. This fact suggests that multiplexity is an important 
ingredient to include in the modelling of real-world systems. 
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The network paradigm has proven to be a successful 
framework to study the intricate patterns of relations 
among the constituents of real-world complex systems, 
from the Internet to the human brain mm, and has re¬ 
vealed that the dynamical behaviours observed in such 
systems, such as information spreading, diffusion, opin¬ 
ion formation and synchronisation, are quite often af¬ 
fected —and to some extent determined— by the struc¬ 
ture of the underlying interaction network |3H6]. How¬ 
ever, the recent availability of massive data sets of so¬ 
cial, technological and biological systems has suggested 
that the classical complex network approach might fall 
somehow short in modelling systems whose elementary 
units can interact through more than one type of connec¬ 
tions. This is typical of many real-world systems, such 
as social networks, where people are connected through 
a variety of social relationships, e.g. kinship, friend¬ 
ship, collaboration, competition, or transportation sys¬ 
tems, which often exploit different communication chan¬ 
nels [MO]. Such systems can be treated in terms of 
multi-layer or multiplex networks EHH], where each 
layer describes a particular type of interaction among the 
nodes of the graph. Some recent works have confirmed 
that multi-layer networks are characterised by new levels 
of complexity m and that the interaction of multiple 
network layers can produce new interesting dynamical 
behaviours p^H2Ql[2^ . 

In the realm of dynamical processes on networks [23] 
the simplicity and -still- the richness of random walks has 
attracted much attention in recent years nils]. Ran¬ 


dom walks are the most simple way to explore a network 
using only local information, and the steady-state prop¬ 
erties of a walk, including characteristic times, limiting 
occupation probability, and coverage, have tight relation¬ 
ships with the structure of the graph upon which the walk 
takes place Eauzi. For this reason, random walks have 
also been successfully used as probes of network proper¬ 
ties, with applications ranging from community detec¬ 
tion [28H3Q] to taxonomy of real-world networks m- 
Moreover, specific flavours of random walks are widely 
used for the exploration of online social networks, infor¬ 
mation networks and the like. 

A class of random walkers of particular interest is that 
of walkers whose motion is biased on the structural prop¬ 
erties of the network [32] . In its simplest possible version, 
the considered biased random walks are Markov processes 
whose transition probability is a parametric function of 
the topological properties of the destination node. In this 
way, by tuning the parameters of the biasing function one 
can force the walk to preferentially visit, or avoid, nodes 
exhibiting high or low values of given topological descrip¬ 
tors, such as the degree, clustering or betweenness. In 
particular, degree-biased random walks have been used 
to define new centrality measures [33l |34| , identify com¬ 
munities [35], and provide optimal exploration of a net¬ 
work using only local information [36]. It has also been 
found that the dynamics of degree-biased random walks 
is strongly affected by the presence of degree-degree cor¬ 
relations in the structure of the network [37H39] , so that 
an appropriate choice of the structural bias can be used 
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to perform efficient sampling of unknown networks. 

In this Article we study several ways in which random 
walks can be extended to multi-layer networks, and we 
show how to devise appropriate ways to bias the walkers 
on the topological properties of the nodes at each layer in 
order to perform an efficient exploration of such systems. 
We notice that random walks have already been applied 
to multi-layer networks, e.g. to quantify the impact of 
failures in interconnected systems [40]. However, we will 
focus here on biased random walks and will investigate 
how the biasing function affects the dispersiveness of the 
walk and the steady-state occupation probability distri¬ 
bution. The aim is to find walks which visit far away 
regions of a multiplex network within a relatively small 
number of steps, a property related to the dispersiveness 
of the walk, and, at the same time, guarantee that the 
probability for a walker to visit any node in the system is 
as close as possible to uniform, thus allowing to sample 
unknown graphs in an efficient way. 

The presence of many interdependent layers allows to 
construct several classes of biased random walks, and 
in particular what we call extensive walks and intensive 
walks, where the difference between the two classes is in 
the dependence of the parameters of the biasing function 
on the number of layers of the system. In the former case, 
the biasing function depends on the structural properties 
of the destination node at all the layers of the system 
(thus, the number of parameters is extensive in the num¬ 
ber of layers), while in the latter case the bias depends on 
intrinsically multiplex properties of the destination node, 
which do not depend explicitly on the number of layers 
of the network. 

For both classes of biasing functions, we provide an¬ 
alytical closed forms for the long-time properties of the 
random walks, in terms of stationary probability distri¬ 
bution and entropy rate lu, and we study the effect of 
different structural properties, including the number of 
layers, the presence and sign of inter-layer degree corre¬ 
lations, the redundancy of edges across layers, the density 
of the multiplex and the heterogeneity of the degree dis¬ 
tributions, on the steady-state behaviour of these walks. 
We find that all these properties have a remarkable effect 
on the maximal dispersiveness and on the steady-state 
occupation probability of biased random walks. 

Finally, we study the diffusion properties of several 
real-world multiplex networks, namely the six continen¬ 
tal airline transportation networks, and we show that in 
those cases the pressure to provide robust route alterna¬ 
tives has somehow hindered the overall diffusion proper¬ 
ties of those systems. 


GENERAL FEATURES OF BIASED RANDOM 
WALKS 

Let us consider a M-layer multiplex network of N 
nodes, i.e. a multi-layer graph in which each node can in¬ 
teract with the other ones by means of M different kinds 
of relationships. A multiplex is fully described by the M- 
dimensional array of the adjacency matrices of its layers 
A = ..., where A^^^^ = G 

and a= 1 if node i and node j are connected at layer a. 
In the following we assume the layers to be unweighted, 
but all the results can be easily extended to to the case 
of weighted multiplexes. 

In general, a random walker on a multiplex is not con¬ 
strained on a single layer and can exploit all the connec¬ 
tions pointing out of the current node, at all layers. A 
synthetic -yet incomplete- description of the topology of 
a multiplex is provided by the overlapping adjacency ma¬ 
trix O = Oij, whose entries Oij = account for the 

total number of connections between two nodes across 
all layers [12]. In particular, we consider the class of 
Markovian random walks defined by the transition prob¬ 
abilities: 


This set up is very general and allows for a variety of dif¬ 
ferent motion rules. In fact, fj can be either a function of 
some topological multiplex properties of the arrival node 
j, or an informative combination of some structural fea¬ 
tures of the destination node, measured at all or at a frac¬ 
tion of the layers. Notice that the unbiased random walk 
on the multiplex is obtained by setting /^ = 1, Vj G U. 
In this case a walker jumps out of node i by travers¬ 
ing one of the edges incident on i chosen with uniform 
probability and independently on the layer to which it 
belongs. It is worth noting that the use of the overlap¬ 
ping adjacency matrix {oij} does not automatically make 
the walk in Eq. § equivalent to a random walk on the 
aggregated graph obtained by flattening all the layers in 
a single network. In general, if the biasing function fj 
depends, either explicitly or implicitly, on the structural 
properties of node j in the multiplex network, the walk 
in Eq. 0 cannot be directly mapped on an equivalent 
walk on the aggregated graph. 

Stationary probability distribution. — Starting from 
the one-step transition probability given in Eq. [^ we de¬ 
rive closed forms for several asymptotic properties of the 
walk. Eollowing an approach similar to that used in 
Ref. [32], we now show that for any choice of the bi¬ 
asing function fj the stationary probability distribution 
P* = {Pi} of biased walks on multiplex networks can 
be analytically derived, under the hypotheses that i) the 
topological overlapping matrix O is primitive and that 
a) fj is a time-invariant function of any property of the 
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destination node j. We start by considering the prob¬ 
ability pi^j{t) that a walker starting at node i will be 
found on node j after exactly t time steps: 

^ X • • • X (2) 

and the dual probability pj^i{t): 

Pj^iit) = ^ X ... X (3) 

Comparing Eq. ([^ with Eq. Q and considering that the 
multiplex is undirected (i.e., Oij = Oji)^ we obtain 

CifiPi^jit) = CjfjPj^iit), Vi, j G E (4) 


where q = matrix O is primitive, 

then a stationary probability distribution exists and 
limt^oo Pi^j {t) = P'j : leading to the expression: 

CifiP* =Cjfjp*. (5) 

By imposing the normalisation condition "^jPj = 1 we 
finally get: 


Cjfi 

'Ei cifi ’ 


(6) 


We notice that Eq. ^ is quite general, since it does not 
explicitly depend on the form of the biasing function or 
on the actual structure of each layer or of the topological 
overlapping matrix O. 

In many real-world application scenarios, e.g. in crawl¬ 
ing the structure of online social networks, it is impor¬ 
tant to guarantee that for long enough times the walk 
will end up visiting all the nodes of the graph with the 
same probability. It is easy to prove that an unbiased 
random walk is not a good choice in this case, since its 
steady-state occupation probability distribution is pro¬ 
portional to the degree sequence, hence an appropriate 
bias should be used to avoid to visit hubs more frequently 
than poorly-connected nodes. In practice, it is not always 
possible to find a walk which produces exactly the same 
stationary occupation probability distribution for all the 
nodes, i.e. p* = p = 1/A^, Vi. However, one could instead 
require that the resulting stationary probability distribu¬ 
tion, although not equal for all nodes, has the minimum 
possible variance. In particular, in the following we will 
focus on the normalised standard deviation of the sta¬ 
tionary probability distribution: 


7]{p*) 


pip*) 


(7) 


where /i(p*) and cr(p*) are the average and the standard 
deviation of p*, respectively. We will look for suitable 
combinations of the parameters of the walk that produce 
the smallest possible value of ?7(p*), corresponding to the 


maximum uniformity of the accessibility of the nodes at¬ 
tainable on a certain multiplex network. 

Entropy rate. — One classical measure to quantify 
the mixedness or dispersiveness of a walk on a graph 
is the entropy rate h = St/t [41], where St is 

the Shannon entropy of the set of all the trajectories of 
length t generated from the walk rule, and h is the min¬ 
imum amount of information necessary to describe the 
process m- In particular, h = 0 only if the walk gen¬ 
erates exactly one possible trajectory, while h is maxi¬ 
mum when all the trajectories are equiprobable. Intu¬ 
itively, walks with a high mixedness can explore remote 
regions of a graph within a relatively small number of 
time-steps. This property is again desirable for the effi¬ 
cient exploration of unknown networks, where only local 
information is available. In particular, it is interesting to 
find a biasing function which guarantees that the walk 
does not remain trapped for too long in any region of the 
graph, and this is usually obtained by maximising the 
dispersiveness of the walk. 

It is possible to show that the entropy rate for a Markov 
process can be expressed as 

InGji), (8) 

which means that h depends only on the walk rule and 
on the stationary probability distribution [32]. By sub¬ 
stituting the analytical expression for p* given in Eq. §) 
into Eq. we get: 


i j i 

(9) 

This expression has a natural upper bound, which reflects 
the case of random walks where all trajectories of the 
same length have equal probability. It is interesting to 
notice that, as shown by Burda et al. in Ref. [42], the 
maximal value of entropy rate attainable by any walk on 
a given single-layer graph depends on the structure of 
the graph, and in particular for an undirected graph it is 
equal to InAmax, where Amax is the maximum eigenvalue 
of the adjacency matrix of the graph. 

This result can be extended to the case of walks on 
multiplex networks as follows. The total number of tra¬ 
jectories of length t generated by a walk defined as in 
Eq. Q is equal to Nf = Ei where is the t-th 

power of the overlapping adjacency matrix. In the limit 
of large t, we have 


Ei Ci/i 


lim 

t^oo 


InNt 

t 


= lnA 


max 5 


( 10 ) 


where A^ax is now the maximum eigenvalue of the over¬ 
lapping adjacency matrix O (this result is a direct con¬ 
sequence of the application of the power method). In 
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general, the maximal value of the entropy rate attain¬ 
able with a particular motion rule will be smaller than 
or at most equal to hmax- Since obtaining high mixed¬ 
ness is a desirable property of a walk in many real-world 
applications, such as when searching for a given resource 
on a graph, in the following we will look for combina¬ 
tions of the parameters of different motion rules which 
can produce high value^of h, to better approximate the 
corresponding value of hmax allowed by the structure of 
the network. 

Heterogeneous mean-field. — In the particular case in 
which the bias function fi depends only on the (vecto¬ 
rial) degree ki = {k\^\ kf\ ..., kf^^} of node i, where by 
definition degree of node i at layer 

a, the expression for the stationary probability distribu¬ 
tion can be considerably simplified. Let us consider a 
heterogeneous mean-field, in which all the nodes belong¬ 
ing to the same degree class k are structurally indistin¬ 
guishable. Under these assumption, and since fi depends 
only on the degree, then for all the nodes i having the 
same degree ki = k we have fi = f^. = /fc, but also 
Ci = Cfc. = Cfc, and similarly: 

Pk= QfkCk = -^fk'^Okk>fk> (11) 


where C is an appropriate normalisation constant to en¬ 
sure that = 1. Eq. (11) means that all the nodes 

in the same degree class will have the same steady-state 
probability of being visited by the walk. Notice that o^k' 
is the expected number of edges connecting two nodes 
whose multiplex degree is respectively equal to k and to 
k'. If we assume that there are no edge correlations, i.e. 
that the probability of having = 1 does not depend 

on the probability of having a= 1 for all the possible 
/3 ^ then we can write: 


P*k = (12) 

k' a=l 


since the expected number Okk' of edges between a node 
with degree k and a node with degree k^ is actually equal 
to the sum of the expected number of edges connecting 
these two nodes at each of the M layers (we indicate 
by the degree at layer a of a node whose vectorial 
degree is equal to k'). If we additionally assume that 
there are no intra-layer correlations, then: 




(fc'[“l) 


(13) 


where P(fc'^“^) is the degree distribution at layer a. In 
the end we find: 


M 


Pk = -^fk'Y^fk' 


a=l 




(14) 


This expression for p^, is quite general, and in particu¬ 
lar it is valid even in the presence of inter-layer degree- 
correlations [43] . Since the heterogeneous mean-field dis¬ 
cards intra-layer and edge correlations, which usually 
contribute to hinder the dispersiveness of a walk, Eq. (14) 
can be readily plugged into the expression of the entropy 
rate in Eq. ^ to obtain an estimate of the maximum 
value of h attainable with a given biasing function on 
a multiplex network with an assigned multiplex degree 
sequence {ki}. 


CLASSES OF BIASED RANDOM WALKS 

The introduction of a biasing function in the motion 
rule is mainly motivated by the necessity to obtain an ex¬ 
ploration of the graph which is more efficient, i.e., faster 
with respect to the time needed to visit all the nodes, 
or more homogeneous, i.e., avoiding heterogeneities in 
the stationary distribution probability, in order to ex¬ 
plore with the same probability each node of the graph. 
In single layer networks these two aims are in general 
antithetical. For instance, a biasing function which max¬ 
imises the mixing of the walk (corresponding to higher 
values of entropy rate) usually produces a quite hetero¬ 
geneous stationary occupation probability, mainly due to 
the fact that a better mixing is obtained by exploiting the 
central role played by hubs. High values of h are usually 
achieved in a single-layer uncorrelated graph by a degree- 
biased walk TTji ^ kj with 6=1, and in general with a 
bias 6 > 0 in graphs with non-trivial degree-degree cor¬ 
relations [32|. On the other hand, a uniform stationary 
occupation probability is obtained by using nji ^ k^ with 
b = —1 in uncorrelated graphs, and in general by a value 
of 6 < 0 for graphs with degree-degree correlations, which 
corresponds to forcing the walkers to preferentially move 
towards poorly connected nodes [38] . 

The richness of multi-layer networks allows the explo¬ 
ration of more complex biasing functions and, as we will 
show in the following, usually produces quite interesting 
dynamics. The reason of such richness is that the multi¬ 
plex degree of a node Hs a vectorial rather than a scalar 
quantity, a fact that allows to construct several degree- 
based biasing functions. In the following we present two 
particular classes of such biasing functions, which we call 
extensive and intensive biases, respectively. 

Extensive bias funetions. — We call extensive those 
walks whose motion rule depends on a function of the 
degrees of the destination node at each of the M layers. 
A first example is that of additive degree-biased walks, 
defined by transition probabilities of the form: 

M 

TTji oc y](fc|“y“ (15) 

a=l 

where G M is the bias exponent associated to layer a. 
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Another example is that of multiplicative degree-biased 
walks, whose transition probability takes the form: 

M 

TTji oc YlikPf- . (16) 

a=l 


We named these walks “extensive” since the number of 
free parameters in the motion rule, namely the exponents 
ba, increases with the number of layers M. This pecu¬ 
liar property of extensive walks allows for a fine-grained 
setting of the bias in order to avoid nodes whose replicas 
on each of the M layers belong to a specific degree class. 
For instance, in the case of a two-layer multiplex, if we 
set bi > 0 and 62 < 0 then the walkers will preferen¬ 
tially move towards node having, at the same time, high 
degree on layer 1 and low degree on layer 2 . It might 
sometimes be desirable for a walker to have such sophis¬ 
ticated motion rules. An example is that of multiplex 
collaboration networks, in which nodes are scientists and 
layers represent co-authorship patterns in different fields. 
In that case, we might use an appropriately biased mul¬ 
tiplex random walk which prefers to move towards nodes 
having a higher degree in a particular field, whose sta¬ 
tionary probability distribution will represent a measure 
of the relative importance of each author in that field. 

However, having a number of parameters which scales 
with the number of layers is not always a desirable prop¬ 
erty, especially if one wants to tune these parameters in 
order to obtain a walk with certain dynamical proper¬ 
ties (e.g., either in terms of stationary probability or in 
terms of entropy rate). This problem is efficiently solved 
by intensive bias functions. 

Intensive bias functions. — We call intensive those 
multiplex walks whose motion bias depends on one or 
more intrinsically multiplex properties of the destination 
node. In the following we will focus on the intensive walk 
whose transition probability reads: 




(17) 


where Oj = J2a overlapping degree of node j 

and Vj is the multiplex participation coefficient of j, and 
is defined as [ 12 ]: 


Vi = 


M 


M -1 



(18) 


We notice that by considering and we are effec¬ 
tively using information about the distribution of the 
edges of the destination node across the layers. In 
particular, for fixed number of layers M, Oi is propor¬ 
tional to the average of the distribution defined by ki = 
, •. •, while Vi gives information about the 

homogeneity of ki^ with 7 ^^ ^ 1 if ^ ^ k\^^ Va 


(i.e., if node i has roughly the same degree at all layers) 
and 7 ^^ ^ 0 if almost all the edges of node i lie on just 
one layer. 

We notice that when bo > 0 the walkers will preferen¬ 
tially move towards hubs, while for 6 ^ < 0 they tend to 
visit the poorly connected nodes more often. Similarly, 
for positive values of bp the walkers will preferentially 
move towards truly multiplex nodes^ i.e. nodes whose 
distribution of edges across the M layers is more homo¬ 
geneous, while for 6 ^ < 0 the walkers prefer to move 
towards focused nodes^ i.e. those having the majority of 
their connections in just one or a few of the M layers m- 
In general, by tuning the two parameters bo and bp we 
can obtain a rich variety of different walks. For instance, 
for bo > 0 and 6^ > 0 , the walkers will be attracted by 
truly multiplex hubs (i.e., nodes with many links, almost 
equally distributed across the layers). Conversely, when 
bo > 0 and 6 ^ < 0 focused hubs are visited often and 
multiplex poorly connected nodes are strongly avoided, 
and so forth. The unbiased multiplex walk is recovered 
for bp = bo = 0 . 

The most interesting characteristic of the intensive 
walk defined by Eq. 0 is that the number of free pa¬ 
rameters is fixed and does not scale with the number of 
layers, as instead happens for extensive walks. We will 
show in the following that intensive walks usually per¬ 
form at least as well as extensive walks, e.g. with respect 
to the maximisation of entropy rate or to the minimisa¬ 
tion of heterogeneity in the stationary occupation prob¬ 
ability distribution. 

It is worth noting that in the case of a duplex, i.e. 
when M = 2, even if the number of biasing parameters 
in intensive and extensive walks is the same, their effect 
on the motion of the walkers is different. Differently from 
bi and 62 , intensive biases do not allow to bias the walk¬ 
ers towards nodes with given properties in a particular 
layer but always consider intrinsically multiplex features, 
such as their total number of connections and their het¬ 
erogeneity. 

In order to explore the differences in the dynamical 
properties (i.e., the entropy rate h and the normalised 
standard deviation of the stationary occupation proba¬ 
bility distribution r]{p'^)) of biased multiplex walks, in 
the top panel of Fig. we report the values of h obtained 
by additive, multiplicative and intensive random walks as 
a function of the two bias exponents in a two-layer mul¬ 
tiplex network whose layers have the same average de¬ 
gree {k) and power-law degree distributions P{k) ^ k~^ 
with 7 = 2.5, with no inter-layer correlations and no 
edge overlap [45|. We notice that also in this simple case 
the three walks have remarkably different behaviours. In 
particular, the additive walk exhibits a relative small sen¬ 
sitivity to the values of the biasing exponents, which re¬ 
sults in smaller variations of h. In fact, there is a large 
region of 61 (i.e. 0 < < 2 ) within which the entropy 

rate is almost constant and not very different from the 
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FIG. 1: Heat-maps of the value of entropy rate h of different multiplex biased walks as a function of the parameters of the 
biasing function. The panels correspond, respectively, to additive [right, (a) and (d)], multiplicative [middle, (b) and (e)] and 
intensive walks [left, (c) and (f)] on uncorrelated duplex networks (in the top panels the two layers have the same power-law 
degree distribution P{k) ~ k~^ with 7 = 2.5, while in the bottom panels the two layers have power-law degree distributions 
with different exponents, namely 71 = 2.2 and 72 = 2.7. In general, the maximum of h is obtained for positive values of the 
two bias parameters, corresponding to extensive walks which move preferentially towards nodes having high degrees on both 
layers, and to intensive walks whose motion rule is biased towards truly multiplex nodes. 


absolute maximum for a relatively large range of values 
of the other exponent 62 , i-e. —5 < 62 < 2 (the same 
reasoning is valid for 0 < 62 < 2 and — 5 < 61 < 2, due 
to the symmetry of the additive bias function). 

Conversely, the picture is much richer and less triv¬ 
ial in the case of multiplicative and intensive walks, for 
which the maximum of h is obtained for a relatively small 
range of parameters, usually corresponding to positive 
exponents. We obtain slightly different results when we 
consider two layers with different power-law degree distri¬ 
butions P{k^^^) ^ and P{k^‘^^) ^ namely 

with exponents 71 = 2.2 and 72 = 2.7 respectively. In 
this case, the symmetry in the additive and multiplicative 
phase diagrams is broken, and the maximum values of h 
are found by biasing the walk towards nodes with high 
degree on both layers, with a higher biasing exponent on 
the degree of the second layer, which has a more homo¬ 
geneous degree distribution. Also the phase diagram for 
the intensive walk is modified, with the line of maximum 
values becoming thinner. 

Similar considerations hold for the phase diagram of 
7 ^(p*), reported in Fig.[^ In this case, the minimum vari¬ 
ance (yielding a more homogeneous exploration of nodes) 
is obtained for negative values of the two bias exponents. 
Moreover, the phase diagram exhibits quite small varia¬ 


tions in the case of additive walk, while we observe more 
heterogeneity in the case of multiplicative and intensive 
walks. Again, the symmetry of the phase diagrams of the 
extensive walks is broken when pairs of layers with dif¬ 
ferent power-law exponents 71 , 72 are considered, with 
the region 62 > hi showing greater variations than for 
62 < . Qualitatively similar differences can be obtained 

with asymmetric layers with respect to other statistical 
properties, such as density. 

All the results for synthetic networks, both in the cur¬ 
rent and following sections, have been obtained for layers 
with N = 10^ nodes and averaged over 1000 realisations. 

HOW THE STRUCTURE OF A MULTIPLEX 
AFFECTS THE WALK 

In this section we illustrate how the structure of the 
multiplex network affects the maximal entropy rate and 
the minimum heterogeneity of the stationary occupation 
probability distribution achievable in the system. 

We focus on five structural aspects, namely i) the pres¬ 
ence and sign of inter-layer degree-degree correlations, ii) 
the existence of edge overlap across layers, in) the num¬ 
ber M of layers of the multiplex, iv) the power-law ex¬ 
ponent 7 of the degree distribution of the layers, and v) 
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FIG. 2 : Heat-maps of the normalised standard deviation of the stationary occupation probability distribution r]{p*) of different 
multiplex biased random walks. Legend as in Fig. In extensive walks, the minimum of rj is always attained for negative 
values of the two exponents, while in intensive walks the minimum of 77 is obtained for bo < 0 and ^ 0 , meaning that walkers 
tend to preferentially move towards nodes with small degrees on both layers. 


their density, measured through the average degree (k). 
Since our focus is on the construction of efficient walks 
(in terms of maximal dispersiveness and of homogeneity 
of the stationary occupation probability) the parameters 
of interest in all the cases are the overall maximum value 
of entropy rate, denoted by /imax, and the minimum value 
of the normalised standard deviation, denoted by r/min, 
obtained by extensive and intensive biased random walks 
as a function of the biasing parameters. 

Effect of inter-layer degree correlations. — In a recent 
work [43] the authors have shown that real-world mul¬ 
tiplex networks are usually characterised by non-trivial 
inter-layer degree correlation patterns. In the same pa¬ 
per the authors propose several methods to quantify the 
presence of inter-layer correlations between a pair of lay¬ 
ers, including the rank correlation among the two degree 
sequences, as measured by the Spearman’s coefficient p. 
If we call the rank of node i due to its degree on 
layer a, the Spearman rank correlation coefficient be¬ 
tween layer a and layer p reads: 



where and are the average ranks of nodes re¬ 
spectively at layer a and layer /3. The coefficient p takes 
values in [—1,1], so that p = 1 if the two degree sequences 


are perfectly correlated (meaning that a hub at layer a 
is also a hub at layer /3), while p = —1 when the two 
degree sequence are perfectly anti-correlated, i.e. when a 
hub on layer a is always a poorly connected node on the 
other layer, and viceversa. Intermediate positive (resp. 
negative) values of p indicate weaker positive (negative) 
inter-layer correlations, while p 0 when the two degree 
sequences are uncorrelated. 

In Fig.|^ a) we report the plot of /imax and Prnin for ex¬ 
tensive and intensive walks on two-layer multiplex net¬ 
works with same average degree and power-law degree 
distributions P{k) ^ k~^ with 7 = 2.5, for different lev¬ 
els of inter-layer degree correlations. As made evident 
by the figure, intensive walks usually perform at least as 
well as extensive walks with respect to both maximisa¬ 
tion of entropy and minimisation of the heterogeneity of 
the stationary occupation probability distribution. This 
suggests that, aside from the actual differences in the 
phase space, intensive walks are able to span the same 
range of values of entropy and 7(79*) by using only two 
parameters, irrespective of the actual numbers of layers 
in the multiplex. 

Effect of edge overlap. — We now investigate the im¬ 
pact of the presence of edge overlap on the long-term dy¬ 
namics of extensive and intensive walks. We recall here 
the definition of overlap for an edge (lj), which is the 
fraction of layers in which the edge (i, j) exists [HI |44], 















3.35 



FIG. 3: Values of /imax (top panels) and ?7min (bottom panels) as a function of the the inter-layer degree correlation coefficient p 
(a), the average edge overlap uj (b) and the number of layers M (c), respectively for additive (triangles), multiplicative (squares) 
and intensive (circles) walks. For the entropy rate, we also show the value of /imax = In A corresponding to the maximum entropy 
random walk (solid line), (a) For all walks, /imax is an increasing function of the inter-layer degree correlation coefficient p, 
and provides a very good approximation of the maximum theoretical entropy rate /imax- Notice that intensive walks perform 
at least as well as the extensive ones, (b) As the overlap increases, the estimates of /imax obtained by the biased walks become 
less precise, while pmin increases as a function of cu. (c) /imax increases and pmin decreases as a function of M. In this case we 
only performed simulations for intensive walks. 


i.e.: 

a=l 

The edge overlap of a multi-layer network is defined as 
the average of ujij over all the node pairs for which Oij 7 ^ 0 
(i.e., for all pairs of nodes which are connected by at least 
one edge): 

.—. ^ 27 2 7 

l,J 

where K is the number of pairs of nodes which are con¬ 
nected in at least one of the M layers. Notice that the 
average edge overlap uj is equal to 1 only if all the M 
layers are identical, while uj = 1/M when every edge is 
present on exactly one of the M layers. 

We started from two-layer multiplex networks obtained 
by coupling identical layers (thus having edge overlap 
equal to 1) with power-law degree distributions P{k) ^ 
k~^ with 7 = 2.5, and then we obtained multiplex net¬ 
works with prescribed values of edge overlap by rewiring 
a certain percentage of the edges of one of the two layers 
in order to maintain the degree sequence unaltered. No¬ 
tice that by construction the resulting multiplex networks 
have maximally positive inter-layer degree correlations 
(i.e., p = 1). As shown in Fig. [^b), ?7niin becomes higher 



FIG. 4: Values of /imax (top panels) and ?7min (bottom panels) 
as a function of the exponent 7 of the the power-law distri¬ 
bution of each layer (a) and of the average degree {k) (b), 
respectively for additive (triangles), multiplicative (squares) 
and intensive (circles) walks. For the entropy rate, we also 
show the value of /imax = In A corresponding to the maximum 
entropy random walk (solid line). As shown, for all walks 
/imax appears to increase as a function of both 7 and {k). 
Smaller variations are also found in the values of ?7min. 

as UJ increases, meaning that higher values of edge over¬ 
lap correspond to a more heterogeneous stationary state 
probability distribution. Conversely, /imax decreases with 
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cj, in accordance with the fact that higher edge overlap 
tends to hinder the dispersiveness of the walk, since a 
smaller number of distinct walks can originate from each 
node. Summing up, multiplex networks having smaller 
values of edge overlap are overall preferable in order to 
maximise the dispersiveness of the walk and to obtain 
a more homogeneous stationary occupation probability. 
In other words, a small edge overlap guarantees a more 
effective exploration of a multiplex network and, at the 
same time, a more homogeneous distribution of the prob¬ 
ability of visiting each node. 

Effect of the number of layers. — It is also interesting 
to study how the dynamical properties of intensive walks 
change when the number of layers M is progressively in¬ 
creased. To this aim, we constructed multiplex networks 
with different number of layers, with no inter-layer degree 
correlations and negligible edge overlap, where all the lay¬ 
ers had power-law degree distributions P{k) ~ k~^ with 
7 = 2.5. As shown in Fig. |^c), hmax is an increasing 
function of M, while r^min decreases as the number of lay¬ 
ers grows. In general, the addition of layers in absence of 
inter-layer correlation flattens the structural differences 
among the nodes of the multiplex, and provides better 
dispersiveness and less heterogeneity in the occupation 
probability distribution. 

Effect of the heterogeneity of the degree distribution. 
— We investigate here how the heterogeneity of the the 
degree distribution of each layer affects hmax and r/min- 
To this aim, we considered pairs of uncorrelated layers 
with the same power-law degree distribution P{k) ^ k~^ 
for different values of 7 , maintaining fixed the average 
degree of the networks (k). The plots in Fig. |^a) con¬ 
firm that both hmax and 77niin grow as 7 increases, i.e. as 
the degree distribution of the layers becomes more ho¬ 
mogeneous. We notice though that the variation in 
appears to be relatively smaller, especially for multiplica¬ 
tive and intensive walks. This result can be explained 
by considering that dispersiveness is favoured by more 
homogeneous degree distributions. Layers with different 
power-law exponents 71 and 72 have been considered in 
the previous section. 

Effect of layer density. — Finally, we focus on the 
effect of layer density, measured through the average de¬ 
gree of the layers (k). Once again we report here the case 
of uncorrelated layers with power-law exponent 7 = 2.5, 
but similar results have been obtained for other values of 
7 . As shown in Fig. |^b), both hmax and r^min increase 
as a function of (k). Layers with different average de¬ 
grees and {k^‘^^) break the symmetry of the phase 

diagrams for h and 77 qualitatively in a similar way as 
pairing layers with different power-law exponents. 

Summing up, the analysis of the impact of structural 
properties on the values of hmax and r^min attainable on a 
multiplex network confirms that positive inter-layer de¬ 
gree correlations, small edge overlap, large number of 
layers, and more homogeneous layers all concur towards 


allowing biased random walks with nearly-optimal dis¬ 
persiveness and closely-to-homogeneous steady-state vis¬ 
iting probability. In other words, a multiplex network 
with a large number of layers and small edge overlap, 
where nodes have roughly the same number of links at all 
layers, can be explored ways more efficiently than a sim¬ 
ilar multiplex network where nodes have disassortative 
inter-layer correlations and edges are redundant across 
layers. 

In the following section we show that the multiplex 
airline transportation networks of all the six continents 
have evolved towards a structure which provides a good 
trade-off between efficient exploration and robustness. 

APPLICATIONS TO REAL-WORLD AIRLINE 
TRANSPORTATION NETWORKS 

As an application, we study here the dynamical prop¬ 
erties of multiplex biased walks on a set of real-world 
systems, namely the six continental airline transporta¬ 
tion networks. In such systems nodes represent airports, 
edges indicate the existence of a route between two air¬ 
ports and each layer is associated to an airline company, 
i.e. all the edges in a layer represent the routes oper¬ 
ated by the corresponding airline. Such networks have 
been introduced and extensively studied in Ref. [43]. As 
shown in Table |T] all such multiplex networks consist of 
a relatively high number of layers. For this reason, we 
will use intensive walks to compute the maximal entropy 
rate hmax and the minimum value of the standard devi¬ 
ation of the stationary distribution r^min • In Table [I] we 
also report for each multiplex the average number of lay¬ 
ers M X uj where each edge exists, the theoretical upper 
value of entropy rate In A, and the values of hmax and 
77min obtained by optimising intensive walks. 


Multiplex 

N 

M 

M X UJ 

In A 

^max 

77 min 

Africa 

238 

84 

1.57 

3.36 

2.20 

1.36 

Asia 

795 

213 

2.16 

4.96 

3.52 

1.17 

Europe 

594 

174 

1.55 

4.60 

3.76 

1.06 

North America 

1029 

143 

1.56 

4.70 

3.75 

1.35 

Oceania 

261 

37 

1.52 

3.71 

2.39 

2.00 

South America 

300 

58 

1.81 

3.66 

2.59 

1.08 


TABLE I: Structural properties of the six continental airline 
transportation systems. For each multiplex, we report the 
number of nodes A, the number of layers M, the average 
number of layers in which an edge exists Mxuj, the theoretical 
upper value of entropy rate InAmax and the extremal values 
hmax and r/min obtained by optimising intensive walks. 

We notice that the efficiency of a transportation sys¬ 
tem is usually measured in terms of the accessibility of 
the locations it serves. In particular, in an ideal trans¬ 
portation system it should be easy to travel between any 
pair of far-apart regions of the network, mostly irrespec¬ 
tive of where exactly those locations are located. Now, 
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discarding the cost associated to the distance between 
the nodes of an airline transportation network, high ac¬ 
cessibility can be obtained by guaranteeing that a trav¬ 
eller can reach remote locations in the system without 
large effort, in terms of number of interchanges, and that 
all locations can be visited with comparable effort. We 
have seen that in the language of random walks these two 
criteria correspond, respectively, to the maximisation of 
dispersiveness and to the minimisation of the standard 
deviation of the visiting probability. 

Hence, we can ask whether the six continental air 
transportation systems can guarantee a good level of 
navigability, i.e. an optimal trade-off between disper¬ 
siveness and homogeneity of the visiting probability. We 
reckon that a more informative analysis of the efficiency 
of these systems would require more detailed information 
about the actual patterns of trips travelled by passen¬ 
gers, the cost associated to each route, the presence of 
non-Markovian effects (people often come back to their 
original place at the end of a trip), the non-stationarity 
of the system due to seasonality, etc.. However, we ar¬ 
gue that biased random walks can still provide useful, 
yet coarse-grained, information about the overall navi¬ 
gability of those systems. Since we cannot modify the 
degree distributions of each of the layers, or the patterns 
of inter-layer correlations, or the actual number of layers 
in each continental air transportation system, we focus 
here in particular on the effects of edge overlap. 

In the previous section we showed that networks with 
high edge overlap uj achieve lower maximal values of dis¬ 
persiveness of the walk and larger heterogeneity of the 
equilibrium occupation probability distribution. When 
two nodes are connected by more than one edge, indeed, 
from a dynamical point of the view some connections are 
wasted, since redundant links do not allow for new paths 
in the network. However, their redundancy might often 
be important for a transportation system, since it makes 
specific connections more robust to single link failures. It 
is not unrealistic to assume that multi-layer transporta¬ 
tion systems from the real-world have developed by sat¬ 
isfying a trade-off between the necessity to provide, at 
the same time, high diffusivity together with reasonable 
levels of robustness. 

Because of the large heterogeneity in the size and num¬ 
ber of layers of the six continental transportation sys¬ 
tems, it is necessary to introduce some kind of normal¬ 
isation which allow to compare the results observed in 
different systems. In order to test the effect of edge over¬ 
lap on the diffusion properties of real-world systems, for 
each of the six multiplex networks we computed the z- 
score of the average edge overlap: 


where {uj) and a{uj) represent respectively the average 
value and the standard deviation of the overlap com¬ 


puted on an ensemble of suitably randomised multiplex 
networks. In particular, for each continental airline sys¬ 
tem we sampled 1000 multi-layer graphs from the config¬ 
uration model which maintains fixed the degree sequence 
of all the layers and rearranges the links on each layer, 
pairing edge stubs at random. We computed also z{hmax) 
and ^^(T^rnin), he. the z-scores of the maximal entropy rate 
and minimum variance over those 1000 multiplex graphs. 

The results reported in Fig. [^confirm that also in real- 
world systems hmax is negatively correlated with edge 
overlap, in agreement with the results obtained on syn¬ 
thetic networks. Similarly, 77max is positively correlated 
with UJ. Notice that we have z{uj) > 0 in all the six conti¬ 
nents, meaning that the edge overlap of the real-world 
systems is always higher than that of the null-model, 
in agreement with the observation that real-world trans¬ 
portation networks tend to guarantee a certain level of 
robustness to failures. However, the quest for robustness 
has a cost in terms of dispersiveness and accessibility. 
In fact, hmax is consistently smaller than the value ob¬ 
served in the randomised systems {z{hmax) < 0) for all 
continents, and similarly the steady-state probability dis¬ 
tribution is consistently larger than that observed in the 
null model (2:(7^niin) >0) 

It is quite interesting to note that the two multiplex 
networks with smallest overlap and overall better diffu¬ 
sion properties are the continental networks of Oceania 
and Europe, which span the least geographical space. We 
can speculate that in such systems some nodes represent¬ 
ing cities in different countries are connected comparably 
well by different modes of transport, such as trains and 
bus, suited for relatively short distances and not included 
in our analysis. This might potentially explain the rel¬ 
ative low number of redundant edges in those two air¬ 
line transportation systems. Conversely, the necessity to 
provide route redundancy has somehow forced the air¬ 
line transportation networks of Asia, South America and 
North America towards slightly less efficient configura¬ 
tions. 


CONCLUSIONS 

In our work we have explored how to extend biased 
random walks to the case of multiplex networks, show¬ 
ing that the richness of multi-layer systems allows to de¬ 
fine several different classes of walks. In particular we 
studied the general features of the so-called extensive 
walkers (where the node properties, as the degree, are 
considered separately at each of the layers with different 
biasing parameters) and intensive walkers (biased on of 
the product two intrinsically multiplex, namely the over¬ 
lapping degree and the participation coefficient) finding 
closed forms for the stationary occupation probability of 
these walks and for the entropy rate, and provided sim¬ 
plified heterogeneous mean-field expressions for the case 
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FIG. 5: z-score of the average edge overlap lu versus the z-scores of the maximal entropy rate hmax (a) and the minimum 
standard deviation of the stationary distribution rjuiin (b) obtained through intensive walks. In agreement with findings for 
synthetic networks, also in real-world systems z{hma^) is negatively correlated with z(uj) - Pearson’s correlation coefficient 
r = —0.70 whereas z(?7min) is positively correlated - r — 0.30, which increases to r = 0.67 excluding the outliner North 
America -. 


in which the multiplex has no correlations. 
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